{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Two-Fluid Stokes implementation\n",
    "Riccardo Rossi, Guillermo Casas\n",
    "UPC BarcelonaTech, CIMNE\n",
    "\n",
    "\n",
    "the Stokes equations are\n",
    "\n",
    "$\\frac{\\partial\\rho u}{\\partial t}-\\nabla\\cdot\\mathbf{C}\\nabla^{s}\\mathbf{u}+\\nabla p=f$\n",
    "\n",
    "$\\frac{\\partial\\rho}{\\partial t}+\\nabla\\cdot\\left(\\rho\\mathbf{\\mathbf{u}}\\right)=0$\n",
    "\n",
    "we consider the following (multiscale) decomposition\n",
    "\n",
    "$\\mathbf{u}=\\mathbf{u_{h}}+\\mathbf{u}_{s}$\n",
    "\n",
    "$p=p_{h}+p_{s}$\n",
    "\n",
    "where $\\mathbf{u_{\\mathbf{h}}}$ belongs to the finite element space.\n",
    "If we now apply the Galerkin method and test the initial equations\n",
    "by $w$ and $q$, we obtain\n",
    "\n",
    "$\\left(\\mathbf{w},f+\\nabla\\cdot\\mathbf{C}\\nabla^{s}\\mathbf{u}-\\nabla p\\right)+\\left(q,-\\frac{\\partial\\rho}{\\partial t}-\\nabla\\cdot\\left(\\rho\\mathbf{\\mathbf{u}}\\right)\\right)=0$\n",
    "\n",
    "substituting the proposed decomposition\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},f-\\frac{\\partial\\rho(\\mathbf{u_{h}+\\mathbf{u}_{s}})}{\\partial t}+\\nabla\\cdot\\mathbf{C}\\nabla^{s}\\left(\\mathbf{u}_{h}+\\mathbf{u_{s}}\\right)-\\nabla(p_{h}+p_{s})\\right) & +\\\\\n",
    "+\\left(q,-\\frac{\\partial\\rho}{\\partial t}-\\nabla\\cdot\\left(\\mathbf{\\rho u}_{h}\\mathbf{+\\rho u}_{s}\\right)\\right) & =0\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "Integrating by parts the convective term\n",
    "\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},\\nabla\\cdot C\\nabla^{s}\\left(\\mathbf{u}_{h}+\\mathbf{u_{s}}\\right)\\right)=-\\left(\\mathbf{\\nabla w},C\\nabla^{s}\\mathbf{u}_{h}\\right)+\\int_{\\Gamma}\\mathbf{w}\\cdot C\\nabla^{s}\\mathbf{u}_{h}\\mathbf{\\cdot n}+ & \\left(\\mathbf{w},2\\nabla\\cdot\\mu\\nabla^{s}\\mathbf{u_{s}}\\right)\n",
    "\\end{align*}\n",
    "Here the second term will only contribute to the external boundary \n",
    "\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},-\\nabla(p_{h}+p_{s})\\right) & =\\left(\\mathbf{\\nabla\\cdot w},p_{h}+p_{s}\\right)-\\int_{\\Gamma}p_{h}\\mathbf{w}\\mathbf{\\cdot n}\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "where we dropped the boundary term on $p_{s}$\n",
    "\n",
    "Developing the continuity equation term and, again, using integration\n",
    "by parts within each element and omitting null boundary terms we\n",
    "obtain\n",
    "\n",
    "\\begin{align*}\n",
    "\\left(q,-\\frac{\\partial\\rho}{\\partial t}-\\nabla\\cdot\\left(\\mathbf{\\rho u}_{h}\\mathbf{+\\rho u}_{s}\\right)\\right) & =\\left(q,-\\frac{\\partial\\rho}{\\partial t}-\\nabla\\cdot\\left(\\mathbf{\\rho u}_{h}\\right)\\right)+\\left(\\nabla q,\\mathbf{\\rho u}_{s}\\right)\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "Putting together all these expressions results in the following equation\n",
    "\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},f\\right)-\\left(\\mathbf{w},\\frac{\\partial\\rho(\\mathbf{u_{h}+\\mathbf{u}_{s}})}{\\partial t}\\right)\\\\\n",
    "-\\left(\\mathbf{\\nabla w},C\\nabla^{s}\\mathbf{u}_{h}\\right)+\\int_{\\Gamma}\\mathbf{w}\\cdot2\\mu\\nabla^{s}\\mathbf{u}_{h}\\mathbf{\\cdot n}+\\left(\\mathbf{\\nabla\\cdot w},p_{h}\\right)+\\left(\\mathbf{\\nabla\\cdot w},p_{s}\\right)-\\int_{\\Gamma}p_{h}\\mathbf{w}\\mathbf{\\cdot n} & +\\\\\n",
    "\\left(q,-\\frac{\\partial\\rho}{\\partial t}-\\nabla\\cdot\\left(\\mathbf{\\rho u}_{h}\\right)\\right)+\\left(\\nabla q,\\mathbf{\\rho u}_{s}\\right) & =0\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "Grouping terms containing the effect of the subscales together we\n",
    "obtain\n",
    "\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},f\\right)-\\left(\\mathbf{w},\\frac{\\partial(\\rho\\mathbf{u_{h}})}{\\partial t}\\right)-\\left(\\mathbf{\\nabla w},\\mathbf{C}\\nabla\\mathbf{u}_{h}\\right)+\\left(\\mathbf{\\nabla\\cdot w},p_{h}\\right)+\\int_{\\Gamma}\\mathbf{w}\\cdot\\mathbf{C}\\nabla^{s}\\mathbf{u}_{h}\\mathbf{\\cdot n}-\\int_{\\Gamma}p_{h}\\mathbf{w}\\mathbf{\\cdot n} & +\\\\\n",
    "+\\left(\\mathbf{w},-\\frac{\\partial(\\mathbf{\\mathbf{\\rho u}_{s}})}{\\partial t}\\right)+\\left(\\mathbf{w},\\nabla\\cdot\\mathbf{C}\\nabla^{s}\\mathbf{u_{s}}\\right)+\\left(\\mathbf{\\nabla\\cdot w},p_{s}\\right) & +\\\\\n",
    "\\left(q,-\\frac{\\partial\\rho}{\\partial t}-\\nabla\\cdot\\left(\\mathbf{\\rho u}_{h}\\right)\\right)+\\left(\\nabla q,\\mathbf{\\rho u}_{s}\\right) & =0\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "We now need to introduce a model for the subscales (**ASGS**)\n",
    "\n",
    "$\\mathbf{u}_{s}=\\tau_{1}\\left(\\mathbf{f}-\\rho\\frac{\\partial\\mathbf{u}_{\\mathbf{h}}}{\\partial t}-\\nabla p_{h}\\right)$\n",
    "\n",
    "$p_{s}=\\tau_{2}\\left(-\\frac{\\partial\\rho}{\\partial t}-\\nabla\\cdot\\left(\\mathbf{\\rho u}_{h}\\right)\\right)$\n",
    "\n",
    "We will make the approximation $\\frac{\\partial\\mathbf{u}_{s}}{\\partial t}\\approx\\mathbf{0}$\n",
    "(quasi-static subscales) and that the fluid is incompressible with\n",
    "uniform density ($\\frac{\\partial\\rho}{\\partial t}=0$ and $\\rho=const$\n",
    "in space)\n",
    "\n",
    "\n",
    "Let's now focus on the remaining stabilization terms one by one\n",
    "\n",
    "If we use linear elements, all higher-than-one order derivatives of\n",
    "finite element functions within each element can be discarded and\n",
    "we have\n",
    "\n",
    "$\\mathbf{\\left(\\mathbf{w},\\nabla\\cdot\\mathbf{C}\\nabla^{s}\\mathbf{u_{s}}\\right)}=-\\mathbf{\\left(\\mathbf{\\nabla w},\\mathbf{C}\\nabla^{s}\\mathbf{u_{s}}\\right)+\\int_{\\Gamma}...}\\approx0.$\n",
    "\n",
    "## STAB PRESS TERM\n",
    "$\\left(\\mathbf{\\nabla\\cdot w},p_{s}\\right)\\rightarrow-\\tau_{2}\\left(\\nabla\\cdot\\mathbf{w},\\left(\\nabla\\cdot\\left(\\mathbf{\\rho u}_{h}\\right)\\right)\\right)$\n",
    "\n",
    "Pressure Equation Terms\n",
    "\n",
    "$\\left(\\nabla q,\\rho\\mathbf{u_{s}}\\right)\\rightarrow\\left(\\nabla q,\\rho\\tau_{1}\\left(\\mathbf{f}-\\rho\\frac{\\partial\\mathbf{u}_{h}}{\\partial t}-\\nabla p_{h}\\right)\\right)$\n",
    "\n",
    "\n",
    "## SYMBOLS TO BE EMPLOYED\n",
    "Shape functions $N_{I}$ and derivatives $\\nabla N_{I}$stored respectively\n",
    "in a vector $\\mathbf{N}$and a matrix $\\mathbf{DN}$\n",
    "\n",
    "we define the following matrices,\n",
    "\n",
    "- $\\mathbf{P}$ such that $\\mathbf{P_{I}}$ is the pressure of node I\n",
    "- $\\mathbf{V}$ such that $\\mathbf{V_{Ik}}$ is the velocity of node I,component K\n",
    "- $\\mathbf{W}$ such that $\\mathbf{W_{I}}=\\mathbf{w_{I}}$\n",
    "- $\\mathbf{Q}$ such that $\\mathbf{Q_{I}=}q_{I}$\n",
    "- $\\mathbf{Vn}$ such that $\\mathbf{Vn_{Ik}}$ is the veocity of node I, component K (OF THE OLD TIME STEP)\n",
    "- $\\mathbf{Vnn}$ such that $\\mathbf{Vnn_{Ik}}$is the veocity of node I, component K (TWO STEPS IN THE PAST)\n",
    "\n",
    "\n",
    "Values on the gauss points are expressed in terms of such variables\n",
    "as\n",
    "\n",
    "- $\\mathbf{acc_{h}:=}bdf_{0}\\mathbf{V}+bdf_{1}\\mathbf{V}_{n}+bdf_{2}\\mathbf{V}_{nn}$a\n",
    "3*1 matrix\n",
    "- $\\mathbf{v_{h}:=\\left(V\\right)^{T}N}$ a 3*1 matrix\n",
    "- $\\mathbf{w_{h}:=W^{T}N}$ a 3*1 matrix\n",
    "- $p_{h}\\mathbf{:=}\\mathbf{P^{T}N}$ a 1*1 matrix\n",
    "- $\\mathbf{q_{h}:=Q^{T}N}$ a 1*1 matrix\n",
    "- $\\mathbf{grad\\_ph}:=\\mathbf{DN^{T}\\cdot P}$ a 3*1 matrix\n",
    "- $\\mathbf{grad\\_vh:=DN^{T}\\cdot V}$ a 3*3 matrix\n",
    "- $\\mathbf{grad\\_qh}:=\\mathbf{DN^{T}\\cdot Q}$ a 3*1 matrix\n",
    "- $\\mathbf{grad\\_wh:=DN^{T}\\cdot W}$ a 3*3 matrix\n",
    "- $\\mathbf{\\epsilon:=grad\\_sym\\_vh:=}\\sum\\mathbf{\\mathbf{B_{I}V_{I}}}$\n",
    "a 3{*}3 matrix (it is the symmetric gradient IN VOIGT FORM!)\n",
    "- $\\mathbf{conv:=\\left(v_{h}\\cdot grad\\_vh\\right)^{T}}$ a 3*1 vector\n",
    "\n",
    "we also introduce the constitutive matrix $\\mathbf{C}$(symmetric\n",
    "6*6 in 3D) and the matrix $\\mathbf{B}$ (6*12 for a tetra) such\n",
    "that $\\sigma=\\mathbf{B}\\epsilon$\n",
    "\n",
    "\n",
    "## IMPLEMENTATION - GALERKIN PART\n",
    "The term\n",
    "\n",
    "$\\left(\\mathbf{w},f-\\frac{\\partial\\rho u_{h}}{\\partial t}+\\nabla\\cdot\\mathbf{C}\\nabla^{s}\\mathbf{u_{h}}-\\nabla p_{h}\\right)$\n",
    "\n",
    "Is implemented on each gauss point as\n",
    "\n",
    "$\\mathbf{w_{h}^{T}}\\mathbf{f-}\\mathbf{w_{h}^{T}\\rho acc_{h}}-\\mathbf{grad\\_sym\\_wh^{T}C}\\mathbf{grad\\_sym\\_vh}+\\mathbf{div\\_}\\mathbf{w_{h}^{T}ph}$\n",
    "\n",
    "\n",
    "## IMPLEMENTATION-STABILIZATION PART\n",
    "\n",
    "**PRESSURE TERMS** $-\\rho\\tau_{2}\\mathbf{div\\_wh}\\mathbf{div\\_vh}+\\mathbf{grad\\_qh^{T}}\\mathbf{\\rho\\tau}_{1}\\mathbf{\\left(\\mathbf{fh-\\rho acc_{h}-grad\\_ph}\\right)}$\n",
    "\n",
    "\n",
    "## IMPLEMENTATION - ENRICHMENT PART\n",
    "here we put terms with penr $p_{enr}$ tested against $\\mathbf{w}$\n",
    "and $p_{h}$ and terms in $\\mathbf{u_{h},p_{h}}$tested against $q_{enr}$\n",
    "\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{\\nabla\\cdot w},p_{enr}\\right)-\\int_{\\Gamma}p_{enr}\\mathbf{w}\\mathbf{\\cdot n} & +\\\\\n",
    "-\\left(\\nabla q,\\rho\\tau_{1}\\nabla p_{enr}\\right)+\\left(\\nabla q_{enr},\\rho\\tau_{1}\\left(\\mathbf{f}-\\rho\\frac{\\partial\\mathbf{u}_{\\mathbf{h}}}{\\partial t}-\\nabla p_{h}-\\nabla p_{enr}\\right)_{s}\\right) & =0\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "which is implemented in the practice as\n",
    "\n",
    "$\\mathbf{div\\_}\\mathbf{w_{h}^{T}penr}+\\mathbf{grad\\_qh^{t}}\\mathbf{\\rho\\tau}_{1}\\mathbf{\\left(\\mathbf{-grad\\_penr}\\right)}+\\mathbf{grad\\_qenr^{T}}\\rho\\mathbf{\\tau}_{1}\\mathbf{\\left(\\mathbf{fh-\\rho acc_{h}-grad\\_ph-grad\\_penr}\\right)}$\n",
    "\\end{document}\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
